function thb
global L m k g
theta0=pi/10;
m=1;k=80;g=9.8;
L0= 1;
L =L0+m*g/k;
[t,u1]=ode45(@thbfun,[0:0.005:15],[L0 0 theta0 0]);
[y1,x1]=pol2cart(u1(:,3),u1(:,1));
y1 =-y1;
figure
ymax = max(abs(y1));
axis([-1.2 1.2 -1.2*ymax 0.2]);
axis off
title('弹簧摆','fontsize',14)
hold on;
R=0.055;
yy =-L0:0.01:0;
xx=R*sin(yy./L0*30*pi);
[a,r]=cart2pol(xx,yy);
a=a+theta0;
[xx,yy]=pol2cart(a,r);
line([-1 1],[0 0],'color','r','linewidth',2)
ball=line(xx(1),yy(1),'color','r','marker','.',...
    'markersize',70,'erasemode','xor');
ball2=line(xx(1),yy(1),'color',[0.5 0.51 0.6],'linestyle','-',...
    'linewidth',1.3,'erasemode','none');
spring=line(xx,yy,'color','g','linewidth',2,'erasemode','xor');
pause(0.5)
for i=1 : length(t)
    yy=-u1(i,1):0.01 :0;
    xx=R*sin(yy./u1(i,1)*30*pi);
    [a,r]=cart2pol(xx,yy);
    a=a+u1(i,3);
    [xx,yy]=pol2cart(a,r);
    set(ball,'XData',x1(i),'YData',y1(i));
    set(ball2,'XData',x1(i),'YData',y1(i));
    set(spring,'XData',xx,'YData',yy);
    drawnow;
end
% subplot(2,1,1);
% plot(t, L0 - u1(:,1), 'b', 'LineWidth', 1.5);
% xlabel('时间 (s)');
% ylabel('径向位移 r');
% title('弹簧摆径向位移随时间变化');
% subplot(2,1,2);
% plot(t, u1(:,3), 'r', 'LineWidth', 1.5);
% xlabel('时间 (s)');
% ylabel('角度位移 \theta');
% title('弹簧摆角度位移随时间变化');
% sgtitle('弹簧摆位移曲线');
% Y_r = fft(L0 - u1(:,1));
% P2_r = abs(Y_r/length(t));
% P1_r = P2_r(1:length(t)/2+1);
% P1_r(2:end-1) = 2*P1_r(2:end-1);
% f_r = (0:(length(t)/2))/max(t);
% [max_amp_r, idx_r] = max(P1_r);
% freq_max_amp_r = f_r(idx_r);
% Y_theta = fft(u1(:,3));
% P2_theta = abs(Y_theta/length(t));
% P1_theta = P2_theta(1:length(t)/2+1);
% P1_theta(2:end-1) = 2*P1_theta(2:end-1);
% f_theta = (0:(length(t)/2))/max(t);
% [max_amp_theta, idx_theta] = max(P1_theta);
% freq_max_amp_theta = f_theta(idx_theta);